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ABSTRACT 

We study the evolution of the mass function of dark matter halos in the concordance ACDM model at high red- 
shift. We employ overlapping (multiple-realization) numerical simulations to cover a wide range of halo masses, 
10 7 - IO'V'Mq, with redshift coverage beginning at z = 20. The Press-Schechter mass function is significantly 
discrepant from the simulation results at high redshifts. Of the more recently proposed mass functions, our results 
are in best agreement with Warren et al. (2005). The statistics of the simulations - along with good control over 
systematics - allow for fits accurate to the level of 20% at all redshifts. We provide a concise discussion of various 
issues in defining and computing the halo mass function, and how these are addressed in our simulations. 

Subject headings: methods: N-body simulations — cosmology: halo mass function 



1. INTRODUCTION 

Dark matter halos occupy a central place in the paradigm 
of gravitationally-driven structure formation arising from the 
nonlinear evolution of primordial Gaussian density fluctuations. 
Gas condensation, resultant star formation, and eventual galaxy 
formation occurs within halos. Consequently, the halo profile 
and mass function are central ingredients in phenomenological 
models of nonlinear clustering of galaxies. The distribution of 
halo masses - the halo mass function - and its time evolution, 
are also sensitive probes of cosmology. 

The halo mass function at the high-mass end (cluster mass 
scales) is exponentially sensitive to the amplitude of the initial 
density perturbations, the mean matter density parameter, fi m , 
and to the dark energy controlled late-time evolution of the den- 
sity field. The last feature, particularly at low redshifts, z < 2, 
allows cluster observations to constrain the dark energy content, 
J7a, and the equation of state parameter, w (Holder et al. 2001). 

The halo mass function is also of considerable interest at high 
redshift, relating to questions such as predictions for quasar 
abundance and formation sites (Haiman & Loeb 2001), the for- 
mation history of collapsed baryonic halos, and the reion- 
ization history of the Universe (Furlanetto et al. 2005). Re- 
cent results from the Wilkinson Microwave Anisotropy Probe 
(WMAP) (Kogut et al. 2003; Spergel et al. 2003) indicate that 
reionization could have begun at redshifts as high as z ~ 20. 
Much of the work on possible reionization scenarios is based 
on the simple Press-Schechter (PS) mass function (Press & 
Schechter 1974, Bond et al. 1991) the use of which can lead 
to imprecise predictions for the reionization history. 

Simulations play a dual role in characterization of the halo 
mass function. If only a few fixed sets of cosmological pa- 
rameters and a finite dynamic range are required, simulations 
can produce valuable results. In order to investigate a variety 
of cosmologies and different scenarios for physical processes, 
e.g., reionization, it is nevertheless very convenient, if not nec- 
essary, to have accurate analytic fitting relations. Simulations 
can be used to validate these fits over a wide (albeit, discretely 
sampled) range of parameters. 

Various numerical studies of the mass function have been 



carried out over different mass and redshift ranges. The clos- 
est to the present work are Reed et al. (2003) and Springel et 
al. (2005); in comparison to their results, our halo mass range 
goes deeper by three orders of magnitude, with good statistics 
and control of systematics out to z = 20, substantially higher 
than in these papers. (We review results from other work be- 
low.) Essentially, the earlier results are in very good agreement 
with the Sheth-Tormen mass function (Sheth & Tormen 1999), 
at redshifts z < 10. As we show below, various fitting formulae 
given in the literature - most tuned to simulation results at z = 
- can differ substantially in their predictions at high redshifts, 
by as much as a factor of two. Therefore, it is important to carry 
out simulations of sufficient dynamic range and accuracy to test 
these predictions. 

In order to extract the mass function from simulations, dif- 
ferent questions have to be addressed, such as: How is the mass 
function to be defined? When do the first halos form in a sim- 
ulation? When must the simulation be started in order to cap- 
ture these halos? What force and mass resolution is required 
to capture halos of a certain mass at a specific redshift? We 
have derived and tested certain criteria to ensure that our simu- 
lations capture the halos of interest; details will be given else- 
where (Lukic et al. 2006). 

The paper is organized as follows. In Section 2 we discuss 
popular mass function formulae, previous work, and strategies 
for determining the halo mass function at high redshifts. The 
simulations and mass function results are discussed in Sec- 
tion 3. Criteria for mass and force resolution and initial redshift 
needed to span the desired mass and redshift range are given 
here. We present our conclusions in Section 4. 

2. THE MASS FUNCTION 

Over the last three decades different fitting functions for the 
mass function have been suggested. The first analytic model 
for the mass function was developed by Press & Schechter 
(1974). Their theory considers a spherically overdense region 
in an otherwise smooth background density field. The over- 
density evolves as a Friedmann universe with positive curva- 
ture. Initially, the overdensity expands, but at a slower rate 
than the background universe (thus enhancing the density con- 
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trast), until it reaches the 'turnaround' density, after which it 
collapses. Although formally this collapse ends with a singu- 
larity, it is assumed that in reality the overdense region will viri- 
alize. For an Einstein-de Sitter universe, the density of such an 
overdense region at the virialization redshift is p « 180p e (z). At 
this point, the density contrast from the linear theory of pertur- 
bation growth [S(x,z) = D + (z)6(x,0)] is 6 c (z) « 1.686(1 +z). For 
fi m < 1, S c (z) evolves differently (see Lacey & Cole 1993), but 
the dependence on cosmology is weak (see e.g., Jenkins et al. 
2001). Thus, we adopt S c = 1 .686 = S C (Q). 

Following the above reasoning and with the assumption that 
the initial density perturbations are given by a Gaussian random 
field, the PS mass function is given by: 



fps(&)=\ exp 
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where a is the variance of the linear density field, f(a,z) = 
(M / p/,){dn / dlna~ l ), and pf, is the background density. 

Using empirical arguments Sheth & Tormen (1999, hereafter 
ST) proposed an improved fit of the following form: 
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with a = 0.707, and p = 0.3. Sheth et al. (2001) interpreted this 
fit theoretically by extending the PS approach to an ellipsoidal 
collapse model. In this model, the collapse of a region depends 
not only on its initial overdensity, but also on the surrounding 
shear field. The dependence is chosen to recover the Zel'dovich 
approximation (Zel'dovich 1970) in the linear regime. A halo 
is considered virialized when the third axis collapses (see also 
Lee & Shandarin 1997). 

Jenkins et al. (2001, hereafter Jenkins) combine high res- 
olution simulations for different cosmologies spanning a mass 
range of over three orders of magnitude [~ (10 12 - lO 15 )/z _1 M0], 
and including several redshifts between z = 5 and z = 0. They 
find that the following fitting formula works exceptionally well 
(within 20%), independent of the underlying cosmology: 

Una" 1 +0.61 1 38 ! . 



(3) 
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The above formula is very close to the nominal ST fit. 

By performing 16 nested- volume simulations Warren et al. 
(2005, hereafter Warren) obtain significant halo statistics span- 
ning a mass range of five orders of magnitude [~ (10 10 - 
10 )A M@]. Their best fit employs a functional form similar 
to an improved version of ST (Sheth & Tormen 2002): 

1.1982 



fwarren(cr) = 0.7234 (a~ l ^ + 0.2538 exp 
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The discrepancy between PS and the more accurate fits is ev- 
ident in Figure 1 where the redshift evolution of the mass func- 
tion is shown. The redshift dependence in the analytic mass 
functions enters only through a(z) = a(Q)d(z), where d(z) is the 
growth factor normalized such that d(0) = 1 . As the functional 
dependence on a is different in the different fits, this leads to 
substantial variation across the fits as a function of redshift. For 
Z = the Warren fit agrees - especially in the low mass range 
below 10 13 M Q - to better than 5% with the ST fit. At the high 
mass end the difference increases up to 20%. The Jenkins fit 
leads to similar results over the considered mass range. Note 
that at higher redshifts and intermediate mass ranges around 
10 9 Mq the disagreement between the Warren and ST fits in- 
creases up to 40%. 



The determination of mass functions at high redshifts is a 
nontrivial task. High redshift halos have very low masses, plac- 
ing heavy demands on the mass and force resolution needed to 
resolve them. These requirements can be achieved in two ways. 
First, a simulation with a very large number of particles and 
high force resolution can be performed. This is expensive, and 
only a very limited number of such simulations can be carried 
out. Second, since determining the mass function is simply a 
question of statistics, many relatively modest simulations with 
moderate particle loading can be performed: this is the strategy 
we adopt here. As simulations can only be trusted until a red- 
shift at which the largest mode is close to becoming nonlinear, 
multiple overlapping box sizes must be used. 

Springel et al. (2005) have recently followed the evolution 
of 2160 3 particles in a 500/r'Mpc box. The high mass and 
force resolution allow them to study the mass function reliably 
out to a redshift of z = 10, covering a mass range of roughly 
\Q m h~ l M & to l0 lb h~ l MQ. Examples of single small-box simu- 
lations include Jang-Condell & Hernquist (2001) (l/i _1 Mpc box 
with 128 3 particles evolved to z = 10) and Cen et al. (2004) 
(4/r'Mpc box, 512 3 particles, evolved to z = 6). Results in 
both papers are claimed to be consistent with PS but without 
detailed quantification. The simulation of Reed et al. (2003) 
is a compromise between the two extremes: a box size of 
50/i _1 Mpc with 432 3 particles and a concomitant halo mass 
range of roughly 10 w h~ 1 M Q to lO 14 - 5 A _1 Af . Reed et al. find 
good agreement (better than 20%) with the ST fit up to z — 10. 
For higher redshifts they find that the ST fit overpredicts the 
number of halos, at z = 15 up to 50%. At this high redshift, 
however, their results become statistics-limited, the mass reso- 
lution being insufficient to resolve the very small halos. 

In this paper we analyze a suite of 50 N-body simulations 
with varying box sizes between 4/r'Mpc and 126/r'Mpc with 
multiple realizations of all boxes to study the mass function at 
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FIG. 1. — Relative residuals of the PS, Jenkins, and Warren mass function 
fits with respect to ST for five different redshifts. Note that the ranges of the 
axes are different in the different panels. We do not show the Jenkins fit below 
masses of 10"/r'Mo since it is not valid in this mass range. 
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redshifts up to z = 20 and to cover a large mass range between 
10 7 h~ l M Q and 10 15 /i _1 M o even at high redshifts. Significantly, 
at z = 20, gas in halos with a mass scale above ~ 10 7 h~ l M & can 
cool via atomic line cooling (Tegmark et al. 1997). 

3. SIMULATIONS AND MASS FUNCTION RESULTS 

All simulations in this paper are carried out with the particle- 
mesh code MC 2 (Mesh-based Cosmology Code). MC 2 has 
been extensively tested against other cosmological simulation 
codes (Heitmann et al. 2005). The chosen values of cosmolog- 
ical parameters are: 

0(01=1-0, fi m = 0.253, n baryon = 0.048, 
cr 8 =0.9, H = 70 km/s/Mpc, (5) 
as set by the latest cosmic microwave background and large 
scale structure observations (MacTavish et al. 2005). The mass 
transfer functions are generated with CMBFAST (Seljak & Zal- 
darriaga 1996). We summarize the different runs, including 
their force and mass resolution in Table 1 . 

We identify halos with the friends-of-friends algorithm (FOF), 
based on finding neighbors of particles at a certain distance (see 
e.g., Einasto et al. 1984; Davis et al. 1985). The halo mass 
is defined simply by the sum of particles which are members 
of the halo. (For connections between different definitions of 
halo masses, see White 2001.) Despite several shortcomings 
of the FOF halo finder, e.g., halo-bridging (see, e.g., Gelb & 
Bertschinger 1994, Summers et al. 1995) or statistical biases 
found by Warren et al. (2005), the FOF algorithm itself is well- 
defined and very fast. 

There are two sources of possible biases in determining in- 
dividual halo masses using FOF. First, the halo may be sam- 
pled with an insufficient number of particles (see Warren et al. 
2005). Second, the effective slope of the halo density profile 
close to the virial radius r V! >, at fixed particle number, also in- 
fluences the FOF mass. If the force resolution of the N-body 
code affects the profile, this too, adds a systematic bias. Here 
we record the mass function for the linking length b = 0.2 FOF 
mass including only the correction of Warren et al. (2005). In a 
follow-up paper (Lukic et al. 2006) we will address systematics 
issues in determining halo masses in detail. 

We now discuss criteria found to be very important for demon- 
strating the convergence and robustness of our results. Details 
will be presented in Lukic et al. (2006). The first issue relates 
to the initial redshift of the simulation. Two conditions are im- 
portant: (i) the simulation must begin sufficiently early that the 
initial ZeFdovich displacement is a small enough fraction of the 
mean interparticle separation A p ; on average a particle should 
not move more than ~ A p /3; (ii) the highest redshift where the 
mass function is to be evaluated must be sufficiently removed 
from the redshift of first-crossing z C wss where particles have the 
first chance to form halos. The stringency of these criteria is 
such that the small boxes require very high starting redshifts, 
e.g., the 4/r'Mpc box had an initial redshift z,„ = 500. This is a 
much earlier starting redshift than those used in previous simu- 
lations; the conventional requirement that all modes in the box 
be linear at the initial redshift proves to be much weaker, and 
therefore inadequate, as a convergence criterion. 

A simple test of how well the simulations track the mass 
function formulae is to follow the number of halos in a spec- 
ified mass bin at a given redshift. For this purpose we con- 
vert the mass function fit into a function of z, defining the halo 
growth function as shown in Figure 2. The evolution of three 
mass bins is shown as a function of z along with results from the 



16/i _1 Mpc boxes. The halo growth function is particularly valu- 
able for determining when the halos at a certain mass should 
first form. This is a good test for problems in simulations aim- 
ing to capture halos with a given minimum mass at some red- 
shift. An example of this is insufficient force resolution in the 
base-grids of adaptive-mesh-refinement (AMR) codes. 

Once the number of particles for a simulation and a desired 
mass for the smallest halo are decided, the required box size is 
fixed. The force resolution needed to resolve the smallest halos 
has then to be determined. Our aim here is not to precisely 
measure the halo profile but simply to be certain that the total 
halo mass is correct. As shown in Heitmann et al. (2005) the 
halo mass is a relatively robust quantity and a simple estimate 
of the force resolution is all that is needed. The force resolution 
must be small compared to the comoving halo virial radius r& 
(with the overdensity relative to the critical density, A ~ 200) at 
all redshifts. The resulting inequality can be stated in the form 



1/3 



(6) 



where 6/ is the force resolution and «/, is the number of particles 
per halo. In the simulations performed here we use a ratio of 
one particle per 64 grid cells, which allows halos with roughly 
50 particles to be captured. It has been shown in Heitmann et 
al. (2005) that this ratio does not cause collisional effects and 
leads to consistent results in comparison with other codes. Mass 
function convergence tests with different force resolutions are 
nicely consistent with the above estimate as shown in Lukic 
et al. (2006); time-step criteria and convergence tests are also 
described there. 

The large set of simulations we have carried out allows us to 
study the mass function at redshifts between z = 20 and z = 0. 
The main results are shown in Figure 3, where the simulation 
data for the mass function are shown along with the Warren, 
PS, and ST fits at different redshifts. At all redshifts the Warren 
fit has the best agreement with the simulations with a scatter of 
approximately 20% and is a numerically significant improve- 
ment over ST. Such a close match is quite gratifying given the 
overall dynamic range of the investigation. The PS fit, over the 
mass range considered, is a poor fit at z > 10, deviating by more 
than a factor of two from the numerical results. 




10 15 
Redshift 

FIG. 2. — Halo growth function for three mass bins for the 16/i~'Mpc box. 
The Warren (solid), ST (long-dashed), and PS (short-dashed) fits are compared 
to simulation data with Poisson error bars. Note the quality of the agreement 
with the Warren fit. 
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Table 1 

Summary of the performed runs 



Mesh 


Box Size 


Resolution 




Z final 


Particle Mass 


Smallest Halo 


# of Realizations 


1024 3 


126/i-'Mpc 


120 /r'kpc 


50 





9.94 x 10 y /r'M 


3.98 x io u hr 


'M 


10 


1024 3 


64 /r'Mpc 


62.5 /r'kpc 


80 





1.30 x 10 9 /r'M 


5.2 x 10 10 h~ 


M 


5 


1024 3 


32 /r'Mpc 


31.25 /r'kpc 


150 


5 


1.63 x 10 s /r'M 


6.52 x 10 9 lr 


l M e 


5 


1024 3 


16 /r'Mpc 


15.63 /r'kpc 


200 


5 


2.04 x 10 7 r'M 


8.16 x 10 s 


'M Q 


5 


1024 3 


8 /r'Mpc 


7.81 /r'kpc 


250 


10 


2.55 x 10 6 /r'M 


1.02 x io 8 /r 


'M 


20 


1024 3 


4 /r'Mpc 


3.91 /r'kpc 


500 


10 


3.19 x 10 5 h~ l M & 


1.27 x 10 7 h~ 


l M 


5 



Note. — Mass and force resolutions of the different runs. The smallest halos we consider contain 40 particles. All simulations are run with 256 3 particles. 
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quired, as White (2002) points out, "it may not be sufficient to 
use a simple parametrized form" in constraining cosmological 
parameters with the mass function. 

The error control criteria developed in this work have a nat- 
ural application in high-resolution AMR simulations in the set- 
ting of refinement and error control criteria. Work in this direc- 
tion is in progress. 
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FIG. 3. — The mass function at 5 different redshifts with Poisson error bars. 
The red line is the Warren fit, blue is Press-Schechter, and black is Sheth- 
Tormen. 



4. CONCLUSIONS AND DISCUSSION 

In this paper we have studied the evolution of the mass func- 
tion starting at a redshift of z = 20 and covering a halo mass 
range of 10 7 to lO 15 /r'M . Our results incorporate new halo- 
based N-body error control criteria that are described in more 
detail in Lukic et al. (2006). We find that the Press-Schechter 
mass function deviates significantly from our results. More re- 
cent mass function fits are in better agreement; in particular, 
the recently introduced fitting function of Warren et al. (2005) 
agrees at the 20% level over the entire redshift range. 

The precise agreement of the numerically obtained halo growth 
function as well as the evolution of the mass function with the 
(evolved z = 0) Warren fit demonstrates the remarkable result 
that the evolution of the mass function is completely controlled 
by the linear growth of the variance of the linear density field. 

In order to find a mass function fit relevant to observa- 
tions, several hurdles remain to be overcome, including reach- 
ing agreement on an appropriate definition of halo mass (White 
2001) and improving the precision and accuracy of N-body 
codes beyond the current state of the art (O'Shea et al. 2005; 
Heitmann et al. 2005). Depending on the level of precision re- 
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